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EVALUATION  AND  ADDITIONAL  DOCUMENTATION 
OF  THE  PARABOLIC  MARCHING  CODE  SURFWAKE 

1.  INTRODUCTION 

In  a  recent  report,  Swean  (Ref.  {!])  presented  results  which  were  obtained  from  the 
numerical  simulations  of  the  wake  of  a  high  speed  surface  ship.  Those  numerical  simu¬ 
lations  were  for  conditions  matching  a  series  of  experiments  which  have  been  conducted 
at  the  David  Taylor  Research  Center  (DTRC)  for  the  surface  wake  of  a  self-propelled 
model  of  a  twin-screw  destroyer.  The  initial  conditions  for  the  calculations  were  gener¬ 
ated  from  the  experimental  data  that  had  been  made  available  to  the  Naval  Research 
Laboratory  (NRL).  Those  calculations  were  performed  with  the  NRL  code  TWAKE. 

Recently,  the  Ship  Wake  Consortium  has  been  formed  by  the  Office  of  Naval  Re¬ 
search  (ONR)  combining  the  cooperative  efforts  of  the  University  of  Michigan.  DTRC 
and  NRL.  As  a  part  of  the  research  work  of  the  consortium,  the  code  SURFW.AKE 
has  been  made  available  to  NRL  for  evaluation  and  documentation.  The  SURFW.\KE 
code  has  several  apparent  advantages  compared  to  the  TW.\KE  code  used  by  Swean 
(Ref.  [1]).  In  the  TWAKE  code,  the  available  free  surface  boundary  condition  is  tin* 
“rigid-lid”  approximation.  SURFWAKE  includes  provisions  for  both  the  ’•rigid-lid" 
boundary  condition  and  a  linearized  free  surface  boundary  condition.  SURFW.AKE 
also  includes  the  capability  for  simulations  of  stratified  flow  fields.  Further.  SURF¬ 
WAKE  has  the  potential  to  be  easier  and  more  economical  to  use  than  TW.\KE.  On 
the  basis  of  these  potentially  attractive  characteristics,  it  is  desirable  to  evaluate  rlu' 
capabilities  of  the  computer  program. 

Unfortunately,  the  version  of  SURFWAKE  retained  by  the  Consortium  is  not 
fully  documented  and  for  this  reason  has  not  been  used  within  Navy  programs  \mril 
this  time.  The  deficiencies  in  documentation  have  to  do  with  how  to  properh’  initiate  a 
computation  with  arbitrary,  perhaps  experimentally  derived,  initial  conditions,  and  thi' 
nature  of  the  options  regarding  the  free  surface  boundary  conditions.  Refs.  [2-4j  jirovidi' 
fairly  extensive  documentation  of  the  program  F.ASTW.A.KE.  a  well-known  comimri'r 
program  within  the  Navy  community.  SURFW.\KE  is  a  derivative  of  this  computer 
program  and  is  described  to  an  extent  in  Ref.  [5j.  Thcut  repoit.  howevt'r.  do('s  not 
contain  documentation  of  the  changes  to  the  original  computi'r  program  whicit  woidd 
be  useful  to  the  programmer  or  u.ser. 

The  purpose  of  this  report  is  three-fold;  the  first  is  to  provide  the  adiiitioiiai 
necessary  documentation  for  initializing  a  surface  slnp  'v,ike  cf^n^o<'.^^f i<in.  ''><•  second 
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is  to  evaluate  the  capabilities  of  SURFWAKE  to  simulate  the  turbulent  wake  of  a 
surface  ship,  and  the  third  is  to  expose  the  nature  of  the  available  linearized  free  surface 
boundary  condition  option  and  to  determine  its  effects  in  a  typiczd  computation.  The 
first  of  these  is  rather  easily  disposed  of  and  that  material  is  included  in  Appendix  A. 
That  material  stemds  eJone  from  the  remainder  of  this  report  aind  the  reader  who  is 
interested  only  in  executing  the  computer  program  for  a  given  fluid  initial  condition 
may  proceed  directly  to  Appendix  A.  A  familiarity  with  the  origintil  FASTWAKE  user’s 
manual  (Ref.  [2])  is,  however,  assumed. 

During  the  course  of  pursuing  the  second  and  third  objectives  above,  some  minor 
programming  errors  were  uncovered.  These  errors  and  the  steps  to  correct  them  will 
be  documented  as  appropriate  in  subsequent  sections  of  this  report.  The  evaluation  of 
the  program’s  ability  to  accurately  simulate  turbulent  ship  wakes  is  less  easily  accom¬ 
plished  due  to  the  lack  of  documented  experimental  data.  The  computations  of  Ref.  [1], 
however,  are  at  least  in  qualitive  agreement  with  preliminary  but  unpublished  exper¬ 
imental  data  (Refs.  [6,7])  and  for  that  reason,  this  report  will  discuss  an  attempt  to 
compute  the  same  flow.  Two  different  free  surface  boundary  conditions  were  employed 
dtiring  the  computations;  the  first  was  the  familiar  “rigid-lid”  boundary  condition  used 
by  nearly  all  computations  of  this  type  (and  in  Ref.  [1|),  and  the  second  was  the  lin¬ 
earized  free  surface  boundary  condition  proposed  by  Meng,  et  al.  [5].  Contrasting  the 
results  of  these  separate  computations  and  with  those  of  Swean  [1]  provides  the  means 
to  accomplish  the  second  and  third  objectives. 

The  next  section  will  give  a  brief  overview  of  the  computational  model  used  in 
SURFWAKE.  Included  will  be  a  discussion  of  the  free  surface  boundary  conditi(ms. 
Where  appropriate,  attention  will  be  drawn  to  the  differences  between  this  code  and 
the  TWAKE  code  referenced  above.  The  final  section  will  present  the  results  of  the 
two  calculations  and  contrast  them  amongst  each  other  and  with  those  of  Ref.  [Ij. 

2.  DESCRIPTION  OF  THE  SURFWAKE  COMPUTATIONAL  MODEL 

The  SURFW.\KE  code  is  a  derivative  of  the  F.\STW.\KE  code  developed  by  Meng 
and  Innis  and  documented  in  Ref.  [2].  The  turbxdence  model  which  was  used  in  that 
version  of  F.A,STW.\KE  used  transport  equations  for  the  turbulence  kinetic  energy  aiul 
the  enstrophy  which  is  the  square  of  the  turbulence  vorticity.  The  turl^ulence  model  was 
subsequently  modified  (Refs.  [3.i{)  to  use  transport  'Equations  for  the  turbulence  kinetic 
energy  amd  the  dissipation  function,  e.  The  governing  equation  set  used  in  F.A.STWAKE 
is  similar  to  the  equation  set  used  in  TW.4KE.  The  SURFW.-VKE  code  was  derived  from 
FASTWAKE  by  Meng  (Ref.  [o|)  by  the  addition  of  boundary  conditions  .suital>le  for  a 
free  surface. 
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The  transport  equations  solved  by  SURF  WAKE  are  the  Reynolds- averaged, 
steady,  parabolic  Navier-Stokes  equations.  The  flow  field  variables  axe  separated  into 
a  non-zero  ensemble  mean  part  and  a  zero-mezin  fluctuating  pairt.  The  mean  field  is 
further  written  as  a  sum  of  an  zunbient  field  term  and  the  departure  from  that  field  due 
to  the  wake.  The  spatial  coordinate  system  is  fixed  in  the  body  so  that  the  ambient 
flow  has  a  constant  axial  component  Uo  along  the  i-axis.  The  x-axis  is  positive  upward 
and  the  y-axis  is  defined  by  the  right-hand-rule.  Thus  the  velocity  and  density  fields 
are  : 


U{x,y,z)=  Uo  +  u{x,y,z)  +  u'(x,y,z) 
=  ^o(-?)  +  v(x,y,z)  -f  v'(x,y,z) 
W(x,y,z)=  w(x,y,z)  +  w'(x,y,z) 

g(x,  y,  z)  =  po(-)  +  p(x,  y,  z)  -f  p'(x,  y,  z) 


The  formulation  allows  for  an  ambient  flow  with  a  trzinsverse  component  Vo(x).  With 
this  velocity  and  density  decomposition,  the  ensemble-averaged  parabolic  Navier-Stokes 
equations  as  solved  by  SURFWAKE  are: 


continuity: 


du  dv  dw 
dx  ^  dy^  dz 


=  0, 


axizd  momentum: 


I 


1 


du  duv  duw  /  d^u  d^u\  d{n'v')  d{u'w') 


vertical  momentum: 


d{Uo  +  u)w  dvw  dww  1  dp  ( d^w  d^w\  d  div'ic']  din-'u-') 

dx  dy  dz  pn  dz  \  dy-  dz-  )  po  Oy  Oz 

transverse  momentum: 


d{Uo  +  u)v  dvv  dvw  1  dp  ( d^v  d^v\  t)(r'rM  Oir'ic'] 

-  -  +  — +  —  = - Z  +  +  —  1  -  - ^ - .  (4) 
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The  turbulence  velocity  and  density  correlations  in  the  above  equations  are  com¬ 
puted  from  the  set  of  algebraic  relations: 

q  +  +  Q))  +  CiC2~s,„  (5) 


^ £  '’’'’O' 


dxj 


dxj  po 


where: 


OV  =  -Crt  f 


T)  _  -7-r  -r-r 


5(j  =  —  +  —  p'"',. 

Po  ■'  Po 

_  ^ 

=  ax,  ax.  ’ 


(6) 

(7) 


(S) 

(9) 

:io) 


and  where  V  =  j’Pfc*  and  Q  =  jQkk  axe  the  production  of  turbulence  kinetic  energy  by 
shear  and  buoyancy,  respectively. 

The  equation  set  is  completed  by  equations  for  the  density  perturbation,  p.  the 
turbulence  kinetic  energy,  q  =  j(u'^  +  v'^  +  w'^),  and  the  isotropic  dissipation  function. 
«: 

density: 


a(Lp  -h  u)p  dvp  drvp  _  dpo  O(v'p’)  d(w'p') 

dx  dy  ^  dz  ^  dz  dy  dz 


ai) 


modelled  turbulence  transport: 


d{Ua  +  u)q  a^  _  ,  ^(r  11^ 

dx  dy  dz  dy\  ’’  e  dy)  dz\  "  e  Oz 


+  P  -  —n-'p'  -  f.  I  121 
po 


modelled  dissipation  function: 


a(f7o  +  u)e  dve  dive  _  a  /  7'  ae\  f  ^  '/* 

dx  dy  ^  dz  dy^  *  e  dy)  dz  \  ^  e  Oz 


q-  Of 


C..-P  -  C./-.  il3' 
'7  '7 
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In  this  form,  the  turbulence  model  contains  the  eight  constants  {Ci,  C2,  C,,  C«, 

Ct2t  Cp\^  Cp2}-  From  Ref.  [4j,  the  values  of  these  constants  are  {0.16,  0.23,  0.15. 
0.10,  1.45,  1.80,  0.21,  1.60}.  The  role  of  these  constants  is  quite  similar  to  that  of 
the  constants  used  in  the  turbulence  model  in  TWAKE;  however,  that  model  uses  four 
constants  in  the  modelling  of  the  kinematic  turbulent  stresses  and  the  model  in  SURF- 
WAKE  uses  only  two.  Also,  the  SURFWAKE  model  has  two  density  related  constants 
and  TWAKE  does  not  solve  a  density  equation. 

It  should  be  noted  that  the  ordering  analysis  used  in  developing  the  equation  set  in 
SURFWAKE  has  resulted  in  the  linearization  of  the  Eixial  momentum  equation  through 
the  omission  of  the  udu/dx  term.  Also  note  that  the  cross-shear  velocity  Vo(x)  has 
been  dropped  in  the  above  representation  of  the  equation  set.  The  density  equation, 
of  course,  reduces  due  to  the  divergence-free  condition  on  velocity  for  constant  density 
flows. 

The  SURFWAKE  code  requires  initial  distributions  of  the  three  components  of 
mean  velocity,  the  turbulence  kinetic  energy,  and  the  dissipation  function  in  order  to 
begin  a  computation.  Additionally,  SURFWAKE  requires  an  initial  distribution  of 
the  density.  One  of  the  principal  objectives  of  the  present  work  is  to  compare  the 
results  of  simulations  of  the  wake  flow  of  the  destroyer  model  obtained  with  SURF¬ 
WAKE  with  the  results  obtained  with  TWAKE.  In  order  to  make  the  comparison  of 
results  meaningful,  the  saune  input  initial  distributions  for  the  three  components  of 
mean  velocity  and  turbulence  kinetic  energy  were  used  for  SURFW.4KE  as  were  used 
for  TWAKE  with  the  addition  of  a  constant  density  distribution. 

As  described  in  Ref.  [1],  the  initial  distribution  for  the  dissipation  function  was 
obtained  from  the  turbulence  kinetic  energy.'  using  Eq.  ( 14)  of  Ref.  [1]  and  an  appropriate 
scale  length  1.  An  analysis  wais  made  of  the  u'v'  equation  in  the  two  turbulence  models 
which  suggested  that  the  value  of  used  in  SURFW.4KE  should  be  approximately 
0.55  times  that  use  in  TWAKE.  This  results  in  l,i  =  0.31  ft  rather  than  0.5S  ft  which 
was  used  in  the  calculations  in  Ref.  [1]. 

In  an  effort  to  further  insure  proper  comparability  of  the  results  between  TW.-\KE 
and  SURFWAKE  ,  four  of  the  constants  in  the  turbulence  model  were  modified  .so  as 
to  match  the  constants  in  TW.4KE  as  nearly  as  possible.  Comparison  of  Eq.  ( 13)  witii 
Eq.  (13)  in  Ref.  [1]  shows  that  C,i  and  C,2  have  essentially  identical  roles  in  the  two 
models.  Thus  C^i  and  C,2  were  set  to  1.44  and  1.92  respectively.  .\lso.  comparison  of 
Eqs.  (12  cind  13)  with  Eqs.  ( 12  and  13)  in  Ref.  [li  shows  that  Cq  in  Eq.  ( 12)  corresponds 
to  C\j(JK  in  Eq.  (12)  of  Ref.  [1]  and  that  C,  in  Eq.  (13)  corresponds  to  Cija,  in  Eq.  (13) 
in  Ref.  [1].  Thus  Cq  and  C,  were  set  to  i).06S  and  0.052  respectively.  Clearly  some 
significant  differences  in  the  turbulence  models  still  exist  and  it  was  not  the  intent  of 
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the  present  study  to  make  the  two  turbulence  models  the  same  but  it  was  intended 
that  features  which  could  be  easily  made  comparable  should  be  comparable. 

The  governing  equations  in  SURFWAKE  are  solved  by  an  Euierian  finite-difference 
method.  The  method  is  one  step  explicit  and  is  second  order  (or  better)  accurate  in  the 
x-direction.  The  space  derivatives  are  center  differenced.  The  finite-difference  method 
uses  a  staggered  mesh  system  with  the  axial  velocity  excess  u,  the  turbulence  kinetic 
energy  q,  and  the  dissipation  function  e  defined  on  cell  centers  but  the  transverse 
velocities  v  and  w  and  the  density  p  defined  on  the  cell  faces.  The  calculation  for  the 
values  of  the  dependent  variables  at  the  new  time  step  (or  x  station)  takes  place  only  on 
the  interior  points  of  the  grid.  The  values  of  the  dependent  variables  on  the  boundary 
grid  points  are  determined  from  the  boundary  conditions  cind  the  variable  values  at  the 
neighboring  grid  points.  Because  the  solution  procedure  uses  a  staggered  mesh  grid, 
the  apphcation  of  the  boundary  conditions  for  the  transverse  velocities  and  the  density 
are  different  in  form  on  the  right  and  left  sides  and  on  the  top  and  bottom  surfaces. 
The  primary  boundary  conditions  are  free-slip  rigid  W2ill  at  the  four  sides. 

As  noted  above,  the  SURFWAKE  code  has  options  for  both  the  “ngid-lid”  approx¬ 
imation  for  the  surface  boundary  condition  as  well  as  a  linearized  free  surface  boundary 
condition.  Free  surface  effects  were  first  included  in  the  F.A.STWAKE  code  by  Meng. 
et  ai.  [5].  When  the  free  surface  option  is  activated,  the  various  flow  variables  have  to 
satisfy  free  surface  boundary  conditions.  The  governing  equations  in  the  fluid  domain 
and  the  three-dimensionad  marching  algorithm  remain  unchanged.  Descriptions  of  the 
different  boundary  conditions  and  their  effects  on  the  solution  are  given  below. 

The  top  row  of  cells  in  the  prograim  SURFW.\KE  is  designated  as  the  free  surface. 
Specifically,  the  variables  aae  required  to  satisfy  the  various  free  surface  boundary 
conditions  on  that  top  row.  The  number  of  grid  points  and  y  -  z  coordinates  of  the 
grid  points  are  not  changed.  By  not  regridding  or  tracking  the  exact  free  surface,  the 
program  implies  a  Taylor  series  expansion  of  the  upper  boundary  equations  about  the 
mean  water  surface  (i,e,  c  =  0).  This  modeling  is  valid  for  waves  with  small  amplitudes 
amd  small  slopes  but  not  for  steep  or  near-breaking  waves. 

The  “rigid- lid”  approximation  is  a  free-slip  rigid  wall  condition  eis  used  at  the 
other  three  sides  of  the  computation  domain.  The  r-gradients  of  the  horizontal  flow 
velocities,  the  turbulence  kinetic  energy,  and  the  dissipation  function  are  set  to  zero  at 
the  boundairy  and  the  vertical  velocity  w  is  .set  to  zero. 

The  significant  turbulence  variables  include  the  Reynolds  stresses,  the  turbuh  nt 
kinetic  energy,  and  the  dissipation  function.  In  addition  to  the  btiundary  «'(qiations 
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that  the  variables  have  to  satisfy,  there  are  also  closure  conditions  on  the  pressure- 
strain  correlation,  the  pressure-density  correlation,  and  the  triple-correlation  term  in 
the  turbulent  kinetic  energy  rate  equation.  These  equations  and  closure  conditions 
have  the  effect  of  changing  the  turbulent  length  scales  to  reflect  the  presence  of  the  free 
surfzw:e.  For  a  more  detailed  justification  of  the  effect  of  turbulence  on  the  free  surface 
boundary  conditions,  refer  to  Meng,  et  al.  [5]. 

The  mean  flow  variables  are  also  required  to  satisfy  free  surface  boundary  condi¬ 
tions.  These  are  typically  expressed  as  kinematic  and  dynamic  relationships.  Given  the 
fixed  grid,  the  problem  is  effectively  linearized  with  respect  to  wave  amplitude.  .\lso. 
the  relatively  large  grid  spacing  precludes  the  numerical  modeling  of  capillary  waves 
and  consequently  the  effects  of  surface  tension  cire  not  considered.  In  the  absence  of 
viscosity  this  model  should,  in  principle,  produce  the  usual  irrotational  linear  gravity 
waves. 

The  modified  free  surface  boundary  equations  of  SURFW.A.KE  are  somewhat  dif¬ 
ferent  than  those  described  by  Meng,  et  al.  [5j.  The  kinematic  condition  for  the  free 
surface  elevation. 


TT  n 

=  w  on  z  =  0 

ax 

and  the  equality  between  the  atmospheric  tind  fluid  pressures. 


(U) 


Pf  =  Pa  =  0  on  r  =  0  (15) 

remain  unchamged.  However,  the  dynamic  condition  for  stress  equilibrium  on  the  sur¬ 
face  has  been  corrected  to  include  tl  e  effects  of  gravity.  By  assuming  that  the  flow  on 
the  surface  is  inviscid  and  steady,  but  not  necessarily  irrotational.  B<?rnotilli's  <'quation 
along  a  streamline  can  be  approximated  by 


Unu  =  —ijT]  on  c  =  0  i  16 i 

Combining  Eqs.  (14)  and  (16)  to  eliminate  the  surface  elevation,  rj.  the  free  surfact' 
boundary  condition  for  the  vertical  velocity  component,  ic.  becomes 
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The  boundzuy  conditions  for  the  horizontal  velocity  components  are  approximated  as 


-^  =  0  and  —  =  0  on  2  =  0.  (18) 

az  dz 

SURFWAKE  has  an  option  which  allows  the  user  to  place  psuedo  wave  dampers  at 
one  or  both  of  the  side  boundaries.  This  option  should  be  activated  when  free  surface 
modeling  is  done,  otherwise  reflected  waves  may  contziminate  the  results.  The  user  is 
cautioned,  however,  that  the  exact  behavior  of  the  dampers  hzis  not  been  determined 
2md  example  cases  should  be  run  for  specific  applications.  Descriptions  of  the  theory 
and  implementation  of  the  dampers  are  given  by  Meng  and  Innis  [2]  and  Chan  [8]. 

3.  RESULTS  AND  DISCUSSION 

The  evaluation  of  the  SURFWAKE  code  was  accomplished  by  making  numerical 
simulations  for  the  wake  of  a  ship  model  for  one  set  of  conditions.  series  of  exper¬ 
iments  have  been  conducted  at  DTRC  to  determine  the  three  dimensional  turbulent 
velocity  field  and  the  free  surface  elevations  in  the  wake  region  of  a  twin-screw  destroyer 
model.  A  line  drawing  of  the  model  is  contained  in  Fig.  1  and  Table  1  displays  relevant 
data  on  the  model  dimensions  and  towing  conditions. 


Table  1.  Model  Ship  Parameters 


Model 

High  Speed  Surface  Ship  ' 

DTRC  No, 

5415-1  1 

Length* 

18.8  ft  , 

Beam 

29.9  in  1 

Draft 

9.8  m 

Prop  Dia. 

8.16  in 

Block  Coef. 

0.506 

Tow  Speed 

4.0  knots  1 

Froude  No. 

0.28 

*  (length  between  forward  and  aft  perpendiculars) 


Swean  [1]  used  the  TW'AKE  code  to  simulate  the  wake  for  two  of  the  experimental 
conditions.  The  experimental  conditic.ts  are  described  in  Ref.  [T  and  will  be  briefly 
summarized  here.  Mean  velocity  and  turbulence  mea,surements  witc  taken  using  both 
laser  doppler  velocimetry  (LDV''  and  hot-film  annemometry  iHFl  r('eimi<ittes.  The 
measurement  domain  consisted  of  a  series  of  vertical  planes  at  various  axial  locations  in 
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the  wake.  A  schematic  of  the  model  and  the  coordinate  system  is  shown  in  Fig.  2a.  and  a 
schematic  of  the  typical  experimental  cross-plane  including  the  major  wake  producing 
elements  of  the  model  is  shown  in  Fig.  2b.  Generally,  the  sensors  were  placed  at  a 
specified  depth  (z)  tmd  measurements  were  taken  point  by  point  at  interv'als  of  2.0  in 
in  the  range  -3.0  <  y  <  33.0.  The  measurement  interval  in  the  r  direction  was  also  2.0 
in  from  -1.6  in  (nearest  to  the  free  surface)  to  a  maximum  depth  of  between  -17.0  in 
to  -28.0  in  depending  on  the  distance  aft  of  the  model. 

A  matrix  showing  the  scope  of  the  measurements  obtained  is  given  in  Table  2. 
The  numbers  in  the  table  refer  to  the  approximate  number  of  points  measured  for  a 
particular  test  which,  is  characterized  by 

1)  the  distance  from  the  aft  perpendicular  (.A.P)  to  the  measurement  plane, 

2)  the  tow  speed, 

3)  the  rotation  direction  of  the  propellers  (or  unpropelled),  and 

4)  the  type  of  sensor  used  for  the  measurements. 

The  elements  denoted  by  an  asterisk  have  been  forwarded  to  NRL  via  Refs.  [6-7]. 

Table  2.  Points/data  planes  from  LDV  and  HF  anemometry 


Feet  from  AP 

-1.7 

-1.1 

-0.4 

6.1 

]  22.0 

Speed,  knots 

3 

a 

3 

6 

4 

6 

a 

3 

■1 

6 

1 

4 

6 

a 

Outboard,  LDV 
Outboard,  LDV 

1 

96 

1 

128 

_ _ 1 

111 

1 

1 

00 

o 

* 

279t 

151 

260* 

1 

i 

167 

Outboard,  HF 

IB 

i 

EH 

■ 

73* 

86* 

68 

■ 

Inboard,  LDV 

i - 

L_.  . 

n 

■ 

■ 

175* 

148 

1 

1  , 

■ 

■ 

■ 

■ 

■ 

182t 

Bi 

*  (forwarded  to  NRL) 

t(denotes  plane  where  symmetry  about  vertical  plane  was  tested) 


Two  of  the  data  sets  referenced  in  Table  2  were  used  to  develop  initial  coiulirious  for  rue 
numerical  simulations  using  TW.A.KE.  These  are  the  4.0  knot.>  data  at  the  10  0  ft  plane 
with  both  outboard  (clockwise)  and  inboard  (counter-clockwise)  propeller  rotation. 

Numericcil  simulations  with  the  SURFW.A,KE  code  were  performed  for  conditions 
corresponding  to  the  inboard  propeller  rotation  data  at  the  tow  speed  of  4.0  knots.  The 
same  data  for  tixial  and  transverse  velocities  and  turbulence  kinetic  enersty  which  weri' 
used  for  the  TW.A.KE  simulations  were  used  for  initializing  the  prc'smit  simulations  with 
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SURFWAKE.  As  described  in  Appendix  A,  modifications  were  made  to  SURFWAKE 
to  permit  the  code  to  accept  empirical  input  data  of  this  type. 

The  input  data  used  by  Swean  [1]  were  for  u,  u,  w  and  g  on  a  uniformly  spaced 
grid  with  0.0  <  y  <  3.0  ft  and  —3.0  <  z  <  0.0  ft  and  with  37  equally  spaced  grid 
points  in  each  direction.  The  calculations  were  made  for  the  starboard  half  plane  of 
the  flow  field  and  a  symmetry  plane  boundary  condition  was  used  on  the  left  boundary 
representing  the  centerplane  of  the  flowfield.  In  the  TWAKE  code,  the  variables  are 
placed  on  the  grid  nodes,  and  the  left  column  of  grid  points  represents  the  plane  of 
symmetry.  In  the  SURFWAKE  code,  the  staggered  mesh  scheme  used  for  locating 
the  dependent  variables  requires  an  extra  column  of  points  on  the  left  boundary  to 
implement  a  symmetry  boundary  condition.  The  actual  plane  of  symmetry  is  the 
second  column  of  grid  points  from  the  left  with  the  left-most  column  of  grid  points 
being  a  column  of  image  points. 

The  simulations  with  SURFWAKE  used  the  same  input  data  used  by  Swean  [1] 
for  u,  V,  w  and  y  on  a  uniformly  spaced  grid  with  —0.08333  <  y  <  3.0  ft  and  -3.0  <  r 
<  0.0  ft  and  with  38  equally  spaced  grid  points  horizontally  amd  37  points  vertically. 
In  the  presentation  of  the  results  obtained  with  SURFWAKE,  the  dependent  ^'ariables 
on  the  image  grid  column  are  shown  for  the  the  aixial  velocity,  the  turbulence  kinetic 
energy,  the  transverse  velocities  and  the  principal  Reynolds  stresses.  The  reader  will 
find  that  Swean  [1]  provides  a  very  complete  description  of  the  initial  plane  data  and 
for  the  variable  distributions  in  the  developing  wake  .  The  discussion  which  follows 
will  thus  be  limited  to  describing  the  differences  between  the  results  of  the  present 
simulations  and  the  differences  between  the  present  results  and  those  of  Ref.  [1], 

The  initial  distributions  of  aixial  excess  velocity  and  turbulence  kinetic  energy, 
normalized  by  appropriate  powers  of  the  free  stream  velocity,  are  shown  in  Figs.  3a 
md  3b.  The  initial  distributions  of  the  transverse  velocities,  displayed  in  the  alternate 
form  of  the  nondimensional  swirl  velocity. 

t/  —  ^  /  2  .  2'v  '/2 

kj  =  —  (e  -I-  u;  ) 

c  0 

are  shown  in  Figs.  3c  and  3d.  The  first  of  these  is  for  the  "rigid-hd"  boundary  roudi- 
tion  and  attention  is  called  to  the  top  row  of  points  for  which  the  velocity  vectors  ere 
all  horizontal.  The  second  is  for  the  linearized  free  surface  boundary  condition  and  it 
should  be  noted  that  on  the  top  row,  the  velocity  vectors  are  not  all  horizontal.  Wit  it 
the  linearized  free  surface  boundary  condition,  the  vertical  velocity  is  deternimed  from 
Eq.  (17).  As  was  noted  above,  a  value  was  chosen  for  the  dissipation  '^cale  Icmitli,  '/• 
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which  gave  a  maximum  value  of  the  Reynolds  stresses  equal  to  the  maximum  vcJue 
in  the  initial  plane  of  the  TWAKE  simulation.  The  Reynolds  stresses  (initial  condi¬ 
tions)  are  computed  in  the  program  from  Eq.  (5)  using  the  experimental  mean  velocity 
and  turbulence  fields.  The  distribution  of  the  resultant  of  the  two  principal  Reynolds 
stresses, 


is  shown  in  Fig.  3e  for  the  initial  computational  plane.  With  the  exception  of  the 
additional  column  of  points  on  the  left  boundary  which  is  required  to  implement  the 
symmetry  plane  boundary  condition  and  the  row  of  surface  points  in  the  linearized 
free  surface  option,  the  velocity  input  data  shown  in  these  figures  are  node  for  node 
identical  to  the  input  data  for  TWAKE  in  Ref.  [ij.  Since  the  Reynolds  stresses  are 
computed  from  the  input  data,  and  the  models  for  the  kinematic  turbulent  stresses  are 
different  in  the  two  computer  codes,  the  distribution  for  is  slightly  different  from 
that  given  in  Ref.  [l]. 

In  the  figures  that  follow,  the  downstream  development  of  the  wake  is  shown 
for  both  types  of  free  surface  boundary  conditions.  Distributions  for  the  velocity, 
turbulence  kinetic  energy  and  principal  Reynolds  stress  are  shown^.for  x  =  16.0  ft 
and  30.0  ft.  The  a)  parts  of  the  figures  show  the  distributions  for  the  '‘rigid-lid” 
boimdary  condition  and  the  b)  parts  show  the  distributions  for  the  linearized  free 
surface  boundary  condition.  For  the  sadce  of  brevity,  in  the  discussion  of  these  figures, 
these  boundary  conditions  will  be  refered  to  tis  the  RL  and  LFS  boundary  conditions. 

The  distributions  for  the  axial  velocity  axe  shown  in  Figs.  4  cind  5  for  16.0  ft  and 
30.0  ft.  In  Fig.  4,  very  little  difference  is  seen  in  either  the  drag  wake  of  the  hull  or 
the  overthrust  wake  wake  of  the  propeller.  .At  x  =  30.0  ft.  the  drag  wake  decayed 
slightly  less  under  the  LFS  boundary  condition  and  the  peak  of  the  overthrust  wake 
moved  slightly  downward  and  outboard. 

The  turbulence  kinetic  energy  distribtitions  are  shown  in  Figs.  6  and  7.  The 
differences  at  16.0  ft  are  very  slight.  The  peak  value  is  slightly  higher  for  the  LFS 
boundary  condition.  .At  30.0  ft  the  differences  are  more  distinct.  The  peak  value  is 
about  5%  higher  for  the  LFS  boundary  condition  and  the  location  of  the  peak  value  is 
somewhat  lower  and  more  outboard  than  for  the  RL  boundary  condition. 

The  distributions  of  swirl  velocity  are  shown  in  Figs.  3  and  9.  .At  16.0  ft.  the  swirl 
patterns  in  the  propeller  wake  are  quite  similar  with  the  max'irnum  swirl  velocity. 
being  about  9%  higher  for  the  LFS  bound.ary  condition.  There  is  little  differ('nce  in 


J 
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the  counter-rotating  vortex  above  the  propeller  wzil^e.  For  the  RL  boundary  condition 
there  is  little  development  of  a  transverse  velocity  field  away  from  the  two  principal 
vortices.  For  the  LFS  boundary  condition,  however,  there  is  a  distinct  outward  cross 
flow  developing  directed  away  from  the  main  wake.  At  30.0  ft,  the  values  for  are 
the  same  to  1%.  For  the  LFS  boundary  condition,  the  outward  cross  flow  is  nearly  as 
strong  as  the  swirl  flow  from  the  propeller,  amd  the  right  boundary  of  the  computational 
domain  is  not  sufficiently  far  away.  The  development  of  the  outwaurd  cross  flow  seems 
to  be  affected  by  the  proximity  of  the  right  edge  of  the  computational  domadn.  As 
the  cross  flow  approaches  the  right  boundary,  its  direction  changes  from  horizontal 
to  vertically  upward,  and  w  is  distinctly  non-zero  on  the  right  boundau'y.  For  both 
botmdairy  conditions,  the  counter-rotating  vortex  above  the  propeller  wake  has  mostly 
disappeared. 

The  distributions  of  are  shown  in  Figs.  10  amd  11.  .4t  16.0  ft.  the  distributions 
for  the  two  boundary  conditions  are  neairly  identical,  except  for  the  top  row.  .\long  the 
top  row,  for  the  RL  boundary  condition,  the  stress  vectors  are  directed  inward  toward  a 
point  above  the  propeller  wake,  and,  for  the  LFS  boundary  condition,  the  stress  vectors 
aure  directed  outward  and  are  in  a  direction  opposite  to  the  direction  of  the  vectors  on 
the  row  below.  At  16.0  ft,  the  maxima  of  the  resultant  of  the  two  principal  Reynolds 
stresses,  Tj^,  are  within  2.5%  in  the  two  simulations.  .4t  30.0  ft.  the  peaik  value  of  Tj.  is 
almost  10%  greater  with  the  LFS  boundary  condition.  The  surface  row  of  points  again 
displays  opposite  directions  for  the  stress  vectors,  and  for  the  RL  boundary  condition 
the  stress  resultant  is  greater  with  a  larger  vertical  component  on  the  two  rows  of  points 
below  the  surface. 

Distributions  of  axiad  velocity  and  turbulence  kinetic  energy  in  the  surface  plane  of 
the  simulations  au’e  shown  in  Figs.  12  and  13.  These  figures  show  the  distributions  for 
the  full  plane  of  the  wake.  The  calculations  were  performed  for  the  starboard  half  of  the 
wake  and  data  from  the  starboard  half  plane  were  reflected  to  the  port  half  so  that  the 
wake  for  the  full  plane  could  be  visualized.  The  velocity  distributions  show  only  small 
differences  for  the  two  boundary  conditions.  The  wake  decays  slightly  more  rapidly 
with  the  RL  boundary  condition.  The  distributions  of  the  kinetic  energy  of  turbulence 
show  that  on  the  surface  the  the  decay  of  the  turbulence  energy  is  more  rapid  with  the 
LFS  boundary  condition.  .4t  30.0  ft  the  level  of  turbulence  kinetic  mergy  in  the  surface 
plane  is  about  30%  less  with  the  LFS  boundary  condition  than  with  the  RL  bountlary 
condition. 

The  longitudinal  development  of  the  maximum  and  minimum  values  of  the  a.xial 
velocity  and  the  turbulence  kinetic  energy  are  shown  in  Figs.  14  and  15  for  the  two 
surface  boundary  conditions.  Fig.  14  shows  distributions  of  the  ma.ximum  axial  velocity 
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in  the  wake,  Um^  the  miniminn  axizJ  velocity  in  the  wake,  Umm,  and  the  maximum 

axial  velocity  in  the  surface  plane,  C/,„.  The  minimum  velocity  in  the  surface  plane 

is  identical  to  the  minimum  velocity  in  the  wzJce.  The  maximum  velocity  in  the  wake 

is  in  the  overthrust  wake,  the  minimum  velocity  in  the  wzike  is  in  the  drag  wake  of  I 

the  hull  and  the  maximum  surface  velocity  is  on  the  wake  centerline  where  there  is  an  H 

overthrusted  region  in  the  initial  plane  data. 

Fig.  15  shows  the  distributions  of  maximum  turbulence  kinetic  energy  in  the  wake, 
qm,  and  on  the  surface,  The  peak  in  the  turbulence  is  seen  in  the  previous  figures 
to  be  between  the  propeller  overthrust  wake  and  the  hull  drag  wake  where  the  gradients 
of  the  axial  velocity  are  the  greatest.  The  location  of  the  maximum  of  the  turbulence  I, 

kinetic  energy  is  consistently  below  the  surface  eind  in  fact  slowly  moves  downward  with 
distance  downstreaun. 

A  close  compairison  of  the  peak  velocity  distributions  shows  no  difference  in  the 
maximum  velocities.  The  minimum  velocities  do  show  some  difference  in  that  with  the 
LFS  boxmdary  condition  the  hull  drag  wake  decays  less  rapidly  with  downsream  distance 
than  with  the  RL  boundary  condition.  For  the  turbulence  kinetic  energy  distributions, 
it  may  be  determined  that  the  maximum  turbulence  in  the  wake  decays  less  with  the 
free  surface  boundary  condition,  whereas,  as  noted  earlier,  the  maximum  turbulence  at 
the  surface  decays  distinctly  more  rapidly  with  the  LFS  boundary  condition. 

The  experimental  data  for  the  downstream  stations  have  not  yet  become  a%-ailable. 

Comparison  may  readily  be  made  with  the  results  of  the  TWAKE  simulations  in  Ref.  [1]. 

Figs.  16a-16d  are  the  results  at  x  =  30.0  reported  by  Swean  [1],  and  they  show  very 
close  agreement  with  the  computed  results  in  the  present  simulation  with  the  “rigid- 
lid”  boundary  condition  (see  Figs.  5a.  7a,  9a,  11a).  The  principal  difference  which  is 
seen  between  the  two  sets  of  results  is  in  the  swirl  velocity  at  x  =  30.0  ft.  The  present 
results  showed  that  the  upper  outboard  rotating  vortex  had  es.sentially  disappeared, 
whereas  the  results  in  Ref.  [1]  still  showed  a  distinct  vorte.x  at  that  location.  In  the 
present  simulation  the  maximum  swirl  velocity  at  30.0  ft  was  0.95%  of  C'o.  while  in 
the  TWAKE  simulation  the  decay  of  the  swirl  was  slower  and  the  maximum  swirl  was 
1.35%  of  Uq.  The  more  rapid  decay  of  may  explain  the  more  rapid  disappearance 
of  the  upper  clockwise  vortex. 

Principal  measures  of  the  wake  development  are  given  in  Table  3a  for  x  =  16. 0  ft 
and  in  Table  3b  for  x  =  30.0  ft.  These  tables  show  the  results  for  the  TWAKE 
simulation  in  Ref.  [1],  and  the  present  results  for  the  “rigid-lid"  and  linc'arizcd  free 
surface  boundary  conditions.  The  decline  of  the  axial  velocities  was  about  the  sam<'  in 
the  three  simulations.  The  peak  of  the  proj)eller  overthrust  wake.  Cm-  d('caye<i  at  riu‘ 
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same  rate  in  the  two  present  simulations  and  more  rapidly  than  in  the  TWAKE  simu¬ 
lation.  The  hull  drag  wake,  Umtn,  decayed  less  rapidly  in  the  SURFWAKE  simulation 
and  least  rapidily  with  the  linearized  free  surface  boimdary  condition.  The  maximum 
of  the  turbulence  kinetic  energy,  qm,  (which  was  located  below  the  surface)  decayed 
at  the  same  rate  in  the  three  simulations.  However,  the  maximum  surface  turbulence, 
,  decayed  more  rapidily  with  the  SURFWAKE  simulations  and  at  about  the  same 
rate  for  the  two  boundary  conditions.  The  maximum  of  the  Reynolds  stresses  as  indi¬ 
cated  by  decayed  most  rapidily  in  the  TWAKE  simulation  and  least  rapidly  in  the 
SURFWAKE  simulation  with  the  linearized  free  surface. 

Table  3.  Comparison  of  Numerical  Results 
a).  16.0  ft 


For  the  “rigid-lid”  boundary  condition  the  differences  in  the  results  are  attributable 
•  to  the  differences  in  the  turbulence  models  used  in  the  two  codes.  Overall,  the  compar¬ 

ison  of  the  three  sets  of  results  show  very  good  agreement  for  the  “rigid- lid"  approxi¬ 
mation  of  the  free  surface  boundary  condition.  In  the  present  simulations,  the  ust*  of 
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the  linearized  free  surface  boundary  condition  has  little  effect  on  the  axial  velocities 
2ind  a  small  effect  on  the  turbulence  kinetic  energy.  There  is  a  distinct  effect  on  the 
development  of  the  transverse  velocities.  This  effect  on  the  evolution  of  the  transverse 
velocities  should  be  studied  further. 

Further  Considerations 

A  basic  purpose  of  this  work  was  to  evaluate  the  suitability  of  the  SURFWAKE 
code  for  use  in  simulating  wakes  of  surface  ships.  During  the  course  of  the  work,  several 
additional  areas  of  concern  developed  and  these  concerns  were  addressed  by  additional 
czdculations.  In  performing  the  simulations  which  have  been  presented  above,  the 
calculations  were  made  only  for  the  starboard  half  plane,  since  the  wake  centerline 
forms  what  should  be  a  very  natural  plauie  of  symmetry.  Ref.  [2j  notes  that  any  of 
the  four  sides  of  the  computational  domain  may  be  set  up  as  a  plane  of  symmetry.  In 
following  the  work  of  Swean  [1]  (and  the  physical  location  of  the  e.xperimental  data), 
the  present  simulations  were  naturally  performed  with  a  plame  of  symmetry  on  the  left 
side.  Calculations  were  perfomed  for  a  full  plzme  case  to  verify  the  appropriateness  of 
approximating  the  full  wake  width  by  calculations  for  only  the  starboard  half  plane. 
The  initial  plane  data  were  appropriately  reflected  across  the  centerplane,  thus  forming 
a  set  of  data  for  the  x  =  10.0  ft  initial  plane  which  had  73  columns  of  data  instead  of 
the  38  columns  for  the  half  plane  case  and  which  was  truly  symmetrical  with  respect 
to  the  centerplane. 

The  transverse  velocity  distributions  for  i  =  30.0  ft  are  shown  in  Fig.  17.  The 
distribution  for  the  full  plane  is  shown  in  Fig.  17a.  Along  the  center  column  and 
directly  above  the  center  of  the  propeller  vortices  can  be  seen  some  evidence  of  local 
asymmetry.  By  expanding  the  scale  of  the  plots,  i.e.  by  plotting  each  half  of  the  full 
plane  flowfield  separately,  differences  in  the  two  sides  can  be  seen  more  clearly.  Fig.  17b 
shows  the  port  half  of  the  flowfield  and  Fig.  17c  shows  the  starboard  half.  In  these  latter 
two  figures,  it  is  more  easily  seen  that  some  of  the  velocity  vectors  in  the  centerplane 
have  a  latertJ  component.  It  may  also  be  seen  that  the  columns  to  each  side  of  the 
centerplane  do  not  quite  match  each  other.  Elsewhere  in  the  flow  field.  a.symmetries 
are  not  readily  seen.  By  considering  the  two  halves  of  the  flowfield  separately,  it  was 
also  determined  that  the  peak  swirl  velocity  was  0.995%  Un  for  the  port  lialf  plane 
and  0.964%  Uo  for  the  starboard  half  plane. 

The  reason  for  cisymmetries  developing  in  the  flow  simulation  is  clearly  in  the  nu¬ 
merical  method.  .4.11  of  the  dependent  variables  are  center-differenced.  However,  most 
of  the  flow  variables  are  assigned  locations  at  the  center  of  the  cells  of  the  computa¬ 
tional  grid  while  the  v  and  w  comp<inents  <if  the  velocity  are  assign<’d  locations  on  the 


faces  of  the  grid  cells.  The  apparent  result  is  that  small  differences  may  be  expected 
to  be  fotmd  between  the  port  and  starboard  halves  of  a  full  pleuie  calculation  when  the 
SURFWAKE  code  is  initiailized  as  it  was  for  these  simulations.  However,  the  overall 
accuracy  of  the  simulation  for  either  half  of  the  flow  field  seems  to  be  quite  satisfactory. 
Comparison  of  Fig.  17c  with  Fig.  9a  shows  only  very  small  differences. 

Another  area  of  concern  arose  during  the  course  of  the  free  surface  ctilculations.  In 
Fig.  9b,  for  the  linearized  free  surface  boundary  condition,  it  was  noted  that  the  trans¬ 
verse  velocities  were  not  sufficiently  small  at  the  right  boundary  of  the  computational 
domain.  Clearly,  simulations  with  the  free  surface  boundary  condition  which  extended 
further  downstream  would  be  contaminated  by  the  transverse  velocity  components  not 
vanishing  on  the  right  side.  The  boundau'y  conditions  for  the  right  side  and  for  the 
bottom  require  that  the  dependent  variables  vanish  on  those  boundaries.  Two  reme¬ 
dies  were  readily  available.  The  code  contained  a  provision  for  a  wave  damper  which 
could  be  easily  activated.  A  slight  change  in  the  coding  was  made  to  permit  the  wave 
damper  to  be  applied  only  at  the  side  that  was  not  being  used  as  a  symmetry  plane. 
The  distribution  of  transverse  velocities  at  i  =  30.0  ft  from  a  simulation  using  the 
wave  damper  is  shown  in  Fig.  18a.  The  wave  dtimper  is  set  to  tiffect  the  five  rightmost 
columns  of  grid  points.  The  damping  is  clearly  effective  in  preventing  the  non-zero 
tranverse  velocities  from  reaching  the  right  boundary.  In  fact,  the  damping  causes  an 
excessively  large  change  in  the  transverse  velocities  between  the  fifth  and  sixth  columns 
of  grid  points.  The  amount  of  the  damping  can  be  adjusted,  but  this  option  was  not 
explored. 

A  second  remedy  for  the  velocity  field  interacting  with  the  right  boundary  is  ex¬ 
pansion  of  the  computational  domsun.  The  domain  is  most  easily  expanded  by  adding 
additional  columns  of  points  to  the  right  side.  A  simulation  was  made  with  double  the 
number  of  grid  points  in  the  latersd  direction.  The  tranverse  velocity  distribution  for 
X  =  30.0  ft  from  that  simulation  is  shown  in  Fig.  ISb.  While  the  flowfield  has  not 
reached  the  outer  boundary,  contamination  of  the  right  boundary  would  be  expected 
to  occvir,  if  the  simulation  were  extended  to  60.0  ft.  The  expanded  grid  required  nearly 
a  doubling  of  computer  time  for  the  simulation. 

The  free  surface  boundary  conditions  shown  in  Eqs.  ( 14)  -  ( IS)  are  correct  for  linear 
gravity  waves.  However,  due  to  the  parabolic  nature  of  the  solution  technique,  it  is  not 
clear  that  the  code  is  capable  of  generating  surface  elevations  that  resemble  the  classic 
Kelvin  wave  pattern.  Certainly  the  trcinsverse  wave  system  will  never  appear.  Using 
the  results  from  the  current  simulation,  the  free  surface  elevation  may  be  calculated  two 
ways.  First,  Eq.  ( 16)  states  that  the  wave  amplitude  is  directly  proportional  to  the  a.xiai 
excess  velocity  component.  Plots  of  the  surface  values  of  the  a.xial  velocity  are  shown 
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in  Fig.  12b  for  veJucs  of  x  from  10.0  to  30.0  ft.  and  in  Fig.  19a  for  values  of  x  from  10.0 
to  190.0  ft.  These  figures  are  chziracteristic  of  a  single  wave  crest  on  each  side  of  the 
plane  of  symmetry,  not  diverging  Kelvin  waves.  The  second  method  for  determining 
the  wave  elevation  is  consistent  with  the  method  of  Ref.  [2]  for  determining  the  value 
of  variables  at  the  n  +  1  time  step  given  their  values  amd  their  derivatives  at  the 
time  step.  For  the  wave  elevation,  this  would  be 


or,  using  Eq.  (14), 


=  q" 


+  Ax 


dx 
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(Ax)^ 
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=  77"  +  Ax—  + 


(Ax)^  dw" 


(19a) 


(196) 


The  results  of  this  calculation  are  shown  in  Fig.  19b  for  values  of  x  from  10.0  to 
190.0  ft.  The  wave  damper  was  active  which  yielded  zero  wave  heights  at  the  outer 
boundary.  Again  the  detzul  of  the  transverse  wave  system  is  lacking.  The  manner  in 
which  SURFWAKE  can  be  used  to  give  realistic  wave  profiles  should  be  the  topic  of  a 
future  study. 

Expansion  of  the  computationeil  domain  by  using  a  variably  space  grid  is  an  option 
that  could  be  considered.  Some  experience  with  a  variably  space  grid  has  shown  that 
extra  care  is  needed  in  obtmning  correct  dependent  variable  values  on  the  surface 
and  symmetry  plane  boundaries.  There  was  also  the  intent  to  maintain  a  high  level 
of  similarity  between  the  present  work  tind  that  of  Swean  [1].  Thus,  the  option  of 
emplo)dng  a  variably  spaced  grid  was  not  explored  in  the  present  work. 


4.  SUMMARY 


The  results  which  were  obtmned  in  the  present  simulations  using  the  SURFW.\KE 
code  with  the  “rigid-lid”  approximation  for  the  free  surface  boundary  condition  agreed 
well  with  the  results  which  Swean  [1]  obtained  with  the  TW.A.KE  code.  The  result.s  for 
the  linearized  free  surface  boundary  condition  agreed  well  for  the  a.xinl  velocities  and 
for  the  turbulence  kinetic  energy.  The  differences  in  the  transverse  velocities  were  large 
but  seemed  to  have  only  small  effects  on  the  evolution  of  the  axial  velocity  distributions 
or  the  distributions  of  turbulence  kinetic  energy.  In  the  SUPEW.4KE  code,  the  axial 
momentum  conservation  equation  is  linearized  and  as  a  result  the  code  is  not  appro¬ 
priate  for  use  when  the  wake  drag  or  overthrust  velocities  are  not  small  with  r('S])ect 
to  the  free  stream  velocity.  SURFW.4.KE  floes  have  a  significant  advantage  compaif'd 
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to  TWAKE  with  respect  to  CPU  requirements.  The  computations  in  the  present  work 
required  on  the  order  of  a  tenth  of  the  computer  times  that  are  required  with  TWAKE. 
This  is  because  of  the  explicit  algorithm  amd  the  linearization  used  in  the  SURFWAKE 
code.  While  the  code  requires  a  smaller  mairching  step-size  for  stability  there  is  no 
requirement  for  the  formation  of  the  large  Jacobian  and  subsequent  matrix  inversion. 
Overall,  the  SURFWAKE  code  seems  quite  suitable  for  wake  calculations  when  the 
“rigid-lid”  approximation  for  the  free  surface  boundary  condition  is  used.  However, 
the  linearized  free  surface  boundary  condition  seems  to  give  poor  simulations  for  the 
transverse  velocities.  Further  work  is  needed  with  the  linearized  free  surface  boundary 
condition. 
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APPENDIX  A  -  Input  data  for  SURFWAKE 


For  the  purposes  of  this  Appendix,  the  modified  code  will  be  distinquished  from  the 
original  program  SURFWAKE  by  the  nzune  SURFWAKE .  v2 .  Many  of  the  input  variable 
names  for  the  program  SURFWAKE .  v2  follow  the  definitions  given  by  Meng  and  Innis  [2] 
and  Meng  et  cil.  [5].  Additions  and  changes  to  the  orginal  I/O  variables  are  noted  here. 
Generally,  lower  caise  characters  are  used  to  identify  the  differences  between  SURFWAKE 
and  SURFWAKE. t2.  For  a  more  complete  description  of  the  parts  of  the  program  that 
axe  unchanged,  the  user  is  referred  to  the  above  references. 

The  program  SURFWAKE .  v2  has  options  that  allow  the  user  to  input  part  or  all  of 
the  program  parameters.  The  initiaJ  plane  data  for  the  wake  can  also  be  calculated 
or  read  in.  This  is  useful  if  experimented  data  exist  for  the  hull-form  being  studied. 
Input  is  read  from  the  various  input  files  (tapes)  in  either  NAMELIST,  formatted  READ 
statements,  or  unformatted  READ  statements.  Output  is  written  to  files  (tapes)  as 
described  below.  The  different  I/O  device  numbers,  the  NAMELIST  variables  and  the 
formatted/unformatted  input  variables  are  described  in  the  following  sections.  The 
output  variables  essentially  follow  the  description  given  by  Meng  and  Innis  [2). 

I/O  Device  Numbers 

Device  number  5  cind  device  number  6  are  reserved  for  the  default  input  and 
output  devices  respectively.  The  other  I/O  devices  in  SURFWAKE. v2  are  .set  in  the 
OPEN  commands  listed  below.  Each  commamd  and  device  number  is  given  with  a  brief 
description. 

open  (unit*  8,file»'<o  he  ujer  jpeci^ed' .form* ‘formatted status* 'old' ) 
Device  8  is  for  fixed  grid,  initial  wake  plane  input  read  in  SUBROUTINE 
DATAIN  tmd  SUBROUTINE  SETUPl.  The  data  are  not  interpolated  onto  a 
new  computational  grid.  The  downstream  grid  remains  unchanged  from 
the  one  read  in. 

open  (unit*10 ,f ile* ‘tapelO ’ ,form=‘unf ormatted' , status* ‘unknown’ ) 

Device  10  contains  summary  information  of  output.  The  variables  are 
written  from  SUBROUTINE  PRINT  and  read  in  SUBROUTINE  RTAPE. 

open  (unit* 11 ,f ile* ‘ tape 11 ’ ,form*‘ unformatted ’ , status* ‘unknown’ ) 

Device  11  contains  current  values  from  the  full  blank  COMMON  block. 

The  variables  are  written  from  SUBROUTINE  WTAPE  and  SUBROUTINE 
PRINT  and  read  in  SUBROUTINE  RTAPE. 
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open  (unit*12,file*‘tapel2' ,form*‘unf ormatted’ ,status=‘ unknown’ ) 
Device  12  contaiins  output  of  the  free  surface  elevation,  ETAA.  The 
values  are  written  in  SUBROUTINE  BOUND?. 

open  (unit*14,file»‘tapel4’ ,fonn*‘formatted’ , status®' unknown ’ ) 
Device  14  contains  output  for  surface  contours  written  to  a  file 
from  SUBROUTINE  PRINTSUR. 

open  (unit»15,f ile*‘tapel5’ ,lorm«‘f ormatted’ , status®' unknown ’ ) 
Device  15  contains  output  of  the  maximum  values  of  the  variables 
from  SUBROUTINE  PRINTMAX  and  SUBROUTINE  DATAIN. 

open  (unit»16,f ile*'tapel6’ ,form«'formatted’ .status® 'unknown’ ) 
Device  16  contains  output  of  the  maximum  values  of  the  variables 
on  the  free  surface  from  SUBROUTINE  PRINTSMAX. 


Input  via  NAMELIST  Statements 


NAMELIST  variables  aire  read  in  on  device  5.  They  are  separated  into  four  groups 
with  groui>-names  of  “namerl”,  “group2”,  “groups”,  and  “group4”.  The  variable 
names,  types,  definitions,  tind  default  values  are  shown  in  Table  I.  In  this  table,  the 
vessel  beam  is  denoted  by  the  symbol  B. 


Table  I.l  NAMELIST  /namerl/  Vauriables 

Variable  Tvve  Definition  Default  Value 

rastrt  Integer  Run  control  flag;  set  equal  to  None 

0  for  new  run.  Set  equtd  to  1  to 
restart  the  program  reading  the 
last  COMMON  blocks  from  tape. 

Restart  information,  the  cycle  number 
from  which  the  problem  is  to  be 
continued. 


IRCYCL  Real 


None 


Table  1.2  NAMELIST  /group2/  Variables 


Variable 

Ims. 

Definition 

Default  Value 

Real 

Grid  parameter:  stretch 
coefficient  in  y-direction. 

0.16334881 

B2 

Real 

Grid  parameter:  stretch 
coefficient  in  z-direction. 

0.14207855 

CONDUC 

Real 

Thermal  difFusivity,  normalized 
by  beam,  B. 

0.0 

CZERQ 

Real 

Maximum  value  of  concentration 
at  wake  center. 

1.0 

DIST 

Real 

Integral  length  scale  for 
turbulence,  U,  normalized  by  B. 

0.22 

DT 

Real 

Initial  mar’ching  increment, 
normalized  by  B. 

0.1 

DVDZSH 

Real 

Transverse  shear  gradient, 
normalized  by  Uo/B. 

0.0 

DYMIM 

Real 

Minimum  mesh  size  in  the  y- 
direction,  normalized  by  B. 

0.05 

DZMIH 

Real 

Minimum  mesh  size  in  the  z- 
direction,  normalized  by  B. 

0.05 

FIUM 

Real 

Internal  “Froude  number"  based 
upon  body  speed,  beam,  and 

Brunt- Vaisala  frequency. 

32.0 

FREE 

Integer 

Free  surface  control  flag;  set 
equal  to  1  for  a  free  surface  on 
z  =  0  otherwise  set  equal  to  0. 

0 
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Variable 

FULL 

I 

I 

G 

ICEN 


Tvve  Definition 


Default  Value 


Integer  Run  control  flag:  0 

0  =  Use  default  values  set  in 

SUBROUTIME  DEFINE  for  input  ; 

1  =  Calculate  subsurface  initial 

plane  wake  profiles  as 
described  by  Ref.  [5]; 

2  =  Calculate  surface  initial 

plane  wake  profiles  as 
described  by  Ref.  [5]; 

3  =  Data  input  on  unit  5  for  u,  v,  w, 

and  q.  These  values  will  be  used 
for  the  interpolated  grid  defined  in 
SUBROUTINE  SETUPl.  SUBROUTINE 
SETUPS  will  be  called  instead  of 
SUBROUTINE  SETUPS. 

4  =  Data  input  on  unit  8  for  u,  v,  w, 

and  q.  The  values  are  read  in 
SUBROUTINE  DATAIN.  They  will  be 
used  on  the  {y,z)  grid  points  read 
in  from  device  8  in  SUBROUTINE 
SETUPl.  SUBROUTINE  SETUPS 
will  be  called  instead  of 
SUBROUTINE  SETUPS.  With  this 
option,  the  grid  will  not  be  rescaled. 

The  input  grid  will  be  used  in  all 
the  downstream  computations. 

Real  Gravitational  acceleration.  12 . 3 

normalized  by  Uq  *  Uo/B  . 

Equivalent  to  the  reciproca’  of 
the  square  of  the  Froude  Number 
based  upon  the  beam. 

Integer  Cell  index  I  closest  to  wake  24 

center. 
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Variable 

IDUM(20) 


IMAX 


Type  Definition  Default  Value 

Integer  Control  flag;  IDUM(l)  =  0 

IDUM(1)...IDUM(19)  determines 
which  variables  will  be  output  in 
line  printer  contour  plots  on 
device  6.  For  M  plots  (M  .LT.  19), 

IDUM(1)...IDUM(M),  should  be 
assigned  values  in  accordance  with 
the  BLOCK  PRINT  identifiers  shown 
below.  IDUM(M+l)  must  be  set  equal 
to  0. 

IDUM(20)  controls  the  number 
of  calls  to  the  SUBROUTINE  INICON. 

When  IDUM(20)  equals  1  the  second 
call  to  SUBROUTINE  INICON  will  be 
skipped.  See  Ref.  [2]  for  the  justification 
behind  calling  INICDN  twice. 


BLOCK  PRINT  identifiers  for  IDUM  variables 
notation  aj  used  in  SUBROUTINE  PRINT 
thus  u'v'  elsewhere  appears  here  as  U*V* 


1  Q 

2  U* 

3  V* 

4  W* 

7  U*V* 

8  U*W* 

9  R*W* 

10  R*W* 

11  RHO* 

12  C 

13  U 

14  V 

15  W 

16  PHI 

17  RS 

18  RHO 

19  E 

20  qs 

21  ES 

22  V*W* 

23  q*v* 

24  q*w* 

25  E*V* 

26  E*W* 

27  C*V* 

2S  C*W* 


Integer  Maximum  number  of  zones  in  the  47 

y-direction  inchiding  fioritious 
zones. 
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Variable  Tvve 
ITEST  Integer 


IZONE  Integer 


JCEM 

JMAX 

MAXCYC 

MFLAG 


Integer 

Integer 

Integer 

Integer 


NFREQB  Integer 


NFREQD  Integer 


Definition 

Flag  to  determine  form  of  selective 
editing  in  the  data  print  out  on 
device  6  after  NFREQE  steps: 

0  =  turbulence  and  mean  flow  data 

1  =  mean  flow  data 

2  =  data  in  scratch  arrays  A1-A5, 

PHI,  RHOEQL,  VSHEAR 

Number  of  zones  included  in  the 
wave  damper  along  the  vertical  edges 
of  the  computationsil  domain. 

Cell  index  J  closest  to  wake 
center. 

Maximum  number  of  zones  in  the 
r-direction  including  fictitious 
zones. 

Maximum  number  of  cycles  allowed 
for  computation. 

Control  flag  set  equal  to  0  to  read 
in  probe  positions  otherwise  set 
equaJ  to  1. 

Number  of  longitudinal  time  steps 
between  line  printer  plots. 

(Equivalent  to  number  of  i-ails  to 
SUBROUTINE  BLOCKP.) 


Default  Value 
0 
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2000 
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2000 


Number  of  longitudinal  time 
between  probe  displacement 
calculations. 


;teps 
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Variable 

Im. 

Definition 

HFREQE 

Integer 

Number  of  longitudinal  time  steps 
between  calls  for  output  by 

SUBROUTINE  EDIT.  Recommend  use 
of  a  large  value  of  NFREQE  to  prevent 
EDIT  outputs. 

NFREQF 

Integer 

Number  of  longitudined  time  steps 
between  writes  of  flow  field 
information  on  the  transverse  plane 
to  tapell. 

NFREQL 

Integer 

Number  of  longitudintd  time  steps 
between  calculation  of  power  law 

NFREQP 

Integer 

Number  of  longitudinal  time  steps 
between  flow  field  short  summary 
output. 

qZERQ 

Real 

Maximum  value  of  turbulence  energy 
at  wake  center,  normalized  by 

Uo  *  Uo. 

T 

Real 

Downstream  distance  at  initial  plane, 
normalized  by  B. 

TEXPND 

Real 

Downstream  distance  at  which  rezone 
will  be  performed,  normalized  by  B. 

TIMEMX 

Resd 

Distance  downstrezun  at  which  run  will 
be  terminated  if  MAXCYC  has  not  been 
reached,  normalized  by  B. 

UD 

Real 

Maximum  axial  velocity  excess 
normalized  by  Uq 

VISCOS 

Real 

Molecular  viscosity  coefficient, 
normalized  by  BUq  so  that 

VISCOS  =  l/ReB. 

VISMAX 

Real 

.Maximum  damping  coefficient  in  wave 
damper  along  vertical  boundaries. 

Default  Value 
1000 

15 

10 

1 

0.00366 

10.0 

2000.0 

300.0 

0.082 

0.0000333 

0.0 


Variable 

Type 

DeHniiion 

Default  Value 

VSWIRL 

Real 

Maximum  swirl,  normalized  by  Uq. 

0.0338 

YMAX 

Real 

Maximum  y-dimension  of  the 
computational  grid.  Used  in  Type-2 
or  Type-3  zoning  in  SUBROUTINE 

SETUPl. 

None 

YZERO 

Real 

Y-distance  from  center  of  cell 
(ICEN,  JCEH)  to  axis  system  fixed 
at  wake  center,  normalized  by 

^Vmin  • 

None 

ZMAX 

Real 

Maximum  z-dimension  of  the 
computational  grid.  Used  in  Type-2 
or  Type-3  zoning  in  SUBROUTINE 

SETUPl. 

None 

ZZERO 

Real 

Z-distance  from  center  of  cell 
(ICEN,  JCEN)  to  axis  system  fixed  at 
wake  center,  normeilized  by 

0.0 

Z12 

Real 

Bottom  boundary  of  first  region 
of  constant  density  stratification. 

100.0 

Z23 

Real 

Bottom  boundary  of  second  region 
of  consteint  density  stratification. 

-100.0 

The  last  two  HAMELISTs,  HAMELIST/group3/  and  NAMELIST/group4/,  are  associ¬ 
ated  with  the  generation  of  initial  plane  data  for  subsurface  and  surface  vessels  respec¬ 
tively.  The  vurious  wake  generating  components  are  identified  and  their  contributions 
linetirly  superimposed  on  the  first  downstream  computational  ijritl.  The  rationalization 
for  this  procedure  and  a  brief  description  of  the  v'ariables  in  the  NAMELISTs  arc'  given 
by  Meng,  et  al.  [5].  They  will  not  be  repeated  here. 

Input  via  READ  Statements 

The  following  four  READ  statements  are  executed  in  SUBROUTINE  CARDSNL: 

read  (5,namerl) 
read  (5,group2) 
read  (5, groups) 
read  (5,group4) 
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They  reference  the  NAMELIST  statements  described  in  the  previous  section.  The  proper 
syntax  rules  for  NAMELIST  input  should  be  followed. 

When  the  NAMELIST/group2/  variable  FULL  equals  1,  2,  or  3,  experimental  data 
may  be  read  in  and  seeded  onto  the  program  generated  grid.  The  experimental  quan¬ 
tities  may  be  one  or  more  of  the  following  one-dimensional  variables; 

RHOAMB,  the  initial  ambient  stratified  density  profile  normalized  to 
be  unity  at  the  wake  center,  or 

VSHEAR,  the  shear  velocity  distribution  normalized  by  the  uniform 
stream  velocity  Uo 

or  one  or  more  of  the  two-dimensional  variables: 

C,  the  dye  concentration, 

Q,  the  turbulence  kinetic  energy  normalized  by  Uo  *  Uo. 

RHO,  the  toted  local  fluid  density  normedized  by  the  dimensional 
value  of  the  ambient  density  at  the  wake  center. 

U,  the  streamwise  perturbation  velocity  normalized  by  ['o, 

V,  the  horizontal  velocity  normalized  by  [’o,  or 

W,  the  vertical  velocity  normalized  by  Uo- 

The  first  line  for  the  experimental  data  is  read  in  SUBROUTINE  DATAIN  by  the 
following  READ  statement: 

READ  (5,1000)  TITLE, IFLAG.NDAT.iprn 
1000  FORMAT (AS, 315) 

where 

TITLE  is  a  descriptive  title  for  the  data  or  ’END"  for  the  last  line  of 
the  experimentjd  input.  This  card  must  be  included  even  if 
FULL  =  0. 

IFLAG  is  the  input  variable  identifier  : 

1  =  RHOAMB 

2  =  VSHEAR 

3  =  C 

4  =  q 

5  =  RHO 

6  =  U 

7  =  V 
S  =  W 


HDAT  is  the  number  of  input  data  values. 

iprn  is  a  variable  to  control  the  printing  of  the  input  values  for  the 
IFLAG  variable  on  device  6. 

Set  iprn  =  1  for  printing,  =  0  for  no  printing. 

If  IFLAG  equals  1  or  2  (one-dimensional  input)  then  the  following  READ  statement 
is  executed: 


READ  (5,1005)  (ZI»(I)  ,DATI1I(I)  ,I»1,NDAT) 

1005  FORMAT (SFIO.O) 

Otherwise,  for  IFLAG  equal  to  3  -  8  (two-dimensional  input)  the  READ  statement  is 
READ  (5,1007)  (YIN(I) .ZIN(I) ,DATIN(I) , 1=1 ,NDAT) 

1007  FORMAKeFlO.O) 

where  (YIM.ZIN)  are  the  {y,z)  coordinates  of  the  IFLAG  variable  and  DATIN  is  the 
actual  value  of  the  variable.  The  data  points  must  be  read  in  sequence  of  ZIN  mono- 
tonicly  increasing  for  one  dimensional  arrays.  For  two  dimensional  arrays,  the  order  is 
arbitrary. 

When  the  HAMELIST/group2/  variable  FULL  equals  4,  experimental  data  for  a  fixed 
grid  and  the  fixed  grid  coordinates  are  to  be  read  in  from  device  S.  The  experimental 
data  is  entered  in  SUBROUTINE  DATAIN  and  the  grid  coordinates  are  entered  in  SUBROU¬ 
TINE  SETUPl.  The  first  lines  in  the  file  attached  to  device  8  contain  the  grid  points. 
They  are  read  in  SUBROUTINE  SETUPl  by  the  following  statements; 

read  (8,*)  imax 
do  i  *  1 ,  imaix 

read(8,*)  y(i) 
end  do 

and 

read  (8,*)  jmax 
do  j  *  1,  jmax 

read(8,*)  z(i) 
end  do 
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t 

I 

The  experimental  quantities  must  be  the  two-dimensionzd  variables;  u,  q,  w,  and  v. 
To  read  them  into  the  program,  the  following  statements  are  executed  in  SUBROUTINE 

DATAIH: 

do  i  «  1 ,  imax 

I  ■  do  j  *  1,  jmajc 

k  ■  k  +  1 

read  (8,690)  u(k) ,q(k) ,w(k) , v(k) 

690  format (4al5 .7) 
end  do 

■  end  do 

The  point  (i,j)  =  (1,1)  corresponds  to  the  lower  left  corner  of  the  grid  and  the  point 
(i,  j)  =  (imax,  jmax)  corresponds  to  the  upper  right  corner.  The  data  are  read  in  row 
by  row. 

After  the  initial  plane  variables  described  above  have  been  read  in  from  device  .5 
®  or  device  8,  the  following  READ  statement  is  executed  in  SUBROUTINE  PRINT; 

READ  (5,15)  (TITLE(I),I»1,20) 

15  FQRMAT(20A4) 

where  TITLE  is  the  identifying  title  of  the  run. 

■  If  NFLAG  =  0  ,  then  an  initizd  call  is  made  to  SUBROUTINE  DIPSCAL  to  read  in  experi¬ 
mental  probe  locations.  SUBROUTINE  DIPSCAL  will  compute  the  isopycnic  displacement 
at  those  locations  [2].  The  locations  are  read  in  on  device  5  with  the  following  READ 
statement: 

READ  (5,1)  NPROBE 

I  1  FORMAT  (110) 

DO  5  N  »  1,  NPROBE 

READ  (5,2)  YPROBE(N),  ZPROBE(N) 

2  FORMAT  (2E10.2) 

5  CONTINUE 

,•  See  Appendix  F  of  Ref.  [2]  for  details  when  using  this  option. 
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(a)  model  coordinate  system 


(a)  contours  of  mean  stream  wise  velocity 
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(c)  swirl  velocity  vectors  —  “rigid-lid”  surface  boundary  condition 


Y 

B 


(d)  swirl  velocity  vectors  -  linearized  free  surface  boundary  condition 


Fig.  3  (Cont  d)  -  Distnbution  of  Initial  Profile  Data  -  t  =  10.0  ft 
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(e)  Reynolds  stress  distribution 
Fig.  3  (Cont'd)  —  Distribution  of  Initial  Profile  Data 
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(a)  “rigid-lid”  surface  boundary  condition 


(b)  linearized  free  surface  boundary  condmon 
Fig.  5  -  Contours  of  mean  streamwise  velocity  x  -  30.0  yr 
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(a)  ‘"rigid-iid”  surface  boundary  condition 
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(a)  “rigid-lid”  surface  boundary  condition 
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(b)  linearized  free  surface  boundary  condition 
Fig.  10  —  Reynolds  stress  distribution  r  =  16.0  ft 
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(a)  "rigid-lid”  surface  boundary  condition 
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(b)  linearized  free  surface  boundary  condition 


Fig.  11  -  Reynolds  stress  distnbution  c  =  30.0  Jr 
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X  -  Xq 
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(a)  “rigid-lid”  surface  boundary  condition 


X  -  X, 


3 

(b)  linearized  tree  surface  boundary  condition 
Fig.  12  -  Contours  of  streainwise  velocity  in  plane  of  free  surface 
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X  -  Xo 


B 

(a)  “rigid-lid”  surface  boundary  condition 


0.0  O  5l5  T5 

X  -  X, 


3 

(b)  linearized  free  surface  boundary  condition 
Fig.  13  —  Contours  of  turbulence  kinetic  energy  in  plane  of  free  surface 
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(c)  full  plane  —  starboard  half  only 


Fig,  17  (Cont’dl  -  Swirl  velocity  vectors;  full  plane  simulation  - 


=  30. 


(b>  double  width  channel 


Swirl  velocity  vectors;  with  wave  damper 
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